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Abstract 

An area-preserving map of the unit sphere, consisting of alternating twists and turns, is 
mostly chaotic. A Liouville density on that sphere is specified by means of its expansion 
into spherical harmonics. That expansion initially necessitates only a finite number of 
basis functions. As the dynamical mapping proceeds, it is found that the number of 
non-negligible coefficients increases exponentially with the number of steps. This is to 
be contrasted with the behavior of a Schrodinger wave function which requires, for the 
analogous quantum system, a basis of fixed size. 
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I. INTRODUCTION 



Chaos is commonly associated with nonlinear dynamics, and the elusiveness of quantum 
chaos is sometimes attributed to the fact that Schrodinger's equation is linear. However, 
the different behaviors of classical and quantum systems cannot be explained so simply. 
Any classical Hamiltonian motion can be described by means of a (possibly singular) 
Liouville density, and the evolution of the latter obeys a linear equation, which can be 
written as 

idf/dt = Lf. (1) 

Here, / is a function of the q n and p n , which are independent variables parametrizing 
phase space, and L is the Liouville operator, or Liouvillian, 




This operator is formally "Hermitian" (over a suitable domain of Liouville functions /) so 
that the time evolution of / is a unitary mapping of phase space. Namely, if there is an- 
o te — funCtion , , * h als0 Sa «e S E, ( 2 ), th e sca lar p roduCt ff 9 U **. 
is invariant in time. This is Koopman's theorem [1]. 

The essential difference between the Liouville equation and the Schrodinger equation is 
that, in the generic non-integrable case, the spectrum of the Liouvillian is continuous [2]. 
This gives rise to fundamental differences between the evolution of Liouville densities 
and that of quantum wave functions for bounded systems that have discrete spectra 
in quantum theory. In particular, any quantum state can always be represented, with 
arbitrary accuracy, by a finite number of energy eigenstates. The time evolution of a 
bounded quantum system therefore is multiply periodic, and will sooner or later have 
recurrences [3] . On the other hand, the Liouville equation involves a continuous spectrum, 
and the evolution of a Liouville density cannot be represented by a finite number of terms. 
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Thanks to this property, a Liouville density is able to become more and more convoluted 
with the passage of time, and it may form intricate shapes with exceedingly thin and long 
protuberances, something that a quantum wave function cannot do [4]. 

The purpose of this work is to apply to the Liouville equation some of the mathematical 
techniques that are standard in quantum mechanics. The Liouville density / is expanded 
into a complete set of orthonormal functions, 

f(q,P) =^c n u n (g,p). (3) 

n 

(Note that these u n are not the eigenf unctions of the Liouvillian — the latter are not 
normalizable since its spectrum is continuous.) The Liouville equation (1) then becomes 

% dc m / dt = Y^ L mn c„, (4) 

n 

where L mn is a Hermitian matrix of infinite order. It follows that \ c n\ 2 is constant in 
time. However, even if we start with a small number of nonvanishing c n , it turns out, as 
will be seen below, that as time passes / spreads over more and more u n . A convenient 
quantitative measure of this spread is the entropy-like expression 

S = ~J2 |Cn| 2 log|c n | 2 . (5) 

n 

This S is analogous to the "entropy" used in the study of quantum chaos [5]. The 
intuitive meaning of S is that e s is a rough indication of the number of basis vectors that 
are appreciably involved in the expansion of / in Eq. (|3|). An appropriate name for S 
could be "dimensional entropy" (or D-entropy). 

Since chaotic systems are characterized by the exponential divergence of neighboring 
trajectories, we expect that the number of basis functions needed for representing / (with 
a given level of accuracy) should also increase exponentially with time. In other words, 
we expect S to increase linearly with time for chaotic systems (in the asymptotic limit of 
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large t). On the other hand, for regular systems, whose trajectories diverge linearly, we 
expect an asymptotically logarithmic growth of S. 

In the following sections, we shall investigate a simple dynamical model that has a 
well documented chaotic behavior, and it will be seen that these guesses are qualitatively 
correct. That model involves a discrete time variable, rather than a continuous time, in 
order to make calculations simpler. Instead of the continuous time evolution (|J), there is 
then a unitary transformation of the components c n (explicitly given below). 



Consider a sequence of mappings of the unit sphere x 2 + y 2 + z 2 = 1, in which each 
step consists of a twist by an angle a around the z-axis (namely, every xy plane turns 
by an angle az), followed by a 90° rigid rotation around the y-axis. The result of these 
consecutive twist and turn operations is 



z' = —x cos(a,2) + y sin(az). 
This map is obviously area preserving. It was extensively investigated, both classically 
and in quantum mechanics, by Haake, Kus and Scharf [6] who called it a "kicked top." 
(It is not really like the motion of a rigid top, because of the torsion.) For low values 
of a, most classical orbits are regular (that is, they are quasi-periodic). As a increases, so 
does the fraction of chaotic orbits, until for a = 3 most of the sphere is visited by a single 
chaotic orbit, as may be seen in Fig. 1. That figure also shows the presence of "forbidden" 
areas, corresponding to regular regions located around fixed points of the map [7]. All 
the following calculations refer to the case a = 3, unless stated otherwise. 

The map (|6]) has interesting symmetries. Given any closed orbit, another closed orbit 
can be obtained by changing the signs of both x and z. This symmetry will be called 



II. THE TWIST AND TURN MAP 



x' = z, 




(6) 
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Ry (it is a rotation by 180° around the y-axis). Moreover, for any pair of distinct orbits 
related by R y , it is possible to obtain a third orbit (of twice the length) by means of a 
180° rotation of that pair of orbits around the rr-axis. These symmetries have important 
consequences for the classification of Liouville densities, as will be seen in the next Section. 

A characteristic property of each orbit is its Lyapunov exponent. In the present case, it 
may be defined as follows: consider an infinitesimal circle drawn on the sphere around the 
initial point of a fiducial orbit. As the map proceeds, this infinitesimal circle is deformed 
into an infinitesimal ellipse, having the same area. The ellipse rotates and stretches or 
contracts in an "erratic" way at each step. Let a n be the length of its semi-major axis 
after n steps. The Lyapunov exponent (per step) is defined as 

A = hm lQg(a " /ao) . (7) 

n— >oo ji 

Such a limit indeed exists for a generic chaotic orbit [8, 9]. In the special case of a regular 
orbit, the pulsations of the ellipse are bounded, and therefore A = 0. 

We have computed in this way the Lyapunov exponents for 10 5 orbits with randomly 
chosen initial points. Each orbit was terminated after 10 4 steps (this usually happened 
when the orbit was regular), or, for chaotic orbits, when the major axis of the ellipse 
exceeded 10 16 (because of the inevitable loss of precision in any further computation). 
The details of the calculations are given in Appendix A, and the results are displayed in 
Fig. 2. About 14% of the orbits are regular. For the chaotic ones, we obtained: 

A = 0.346 ± 0.071 (average ± standard deviation). (8) 

Ideally, A should have had a sharp value, the same for all chaotic orbits. This dispersion 
free value could in principle have been found by more sophisticated methods [10], but it 
was not needed for our present purpose. 
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III. MAPPING OF LIOUVILLE DENSITIES 



Instead of considering individual orbits, that is, mapping of points into points, a more 
general approach, which gives new and instructive insights, is the mapping of Liouville 
densities. Let us imagine that an infinitesimal "mass" pdA (which may be positive or 
negative) is attached to each area dA. Let us further assume that this mass is conserved 
during the twists and turns of the unit sphere, so that its density p obeys the linear law 

P'(x',y',z') = p(x,y,z), (9) 

by virtue of dA' = dA. Note that, as the mapping @ is continuous, nearby points are 
mapped into nearby points, so that a blob is always mapped into a single blob (never into 
several disjoint blobs). 

For the twist and turn map ([]), these Liouville densities may belong to three invariant 
symmetry classes, according to their behavior under R y and R x (namely, 180° rotations 
around the y- and x-axes, respectively). For example, if p = F(x 2 , y 2 , z 2 , xyz) is a single 
valued function of its four arguments, this p is even under R y and R x , and is mapped by 
Eq. (9) into another function of the same type. Likewise, any p = yF(x 2 ,y 2 , z 2 ,xyz) is 
even under R y and odd under R x , and is mapped by (9) into a function of the same type. 
For instance, the function p = y has this property, because 

y — > y' = x sin(a,2) + y cos(az) = y 

In general, any p(x, y, z) can be written as the sum of three terms, belonging to one 
the symmetry classes listed in Table 1. 



cos \az) 



xyz sin(az) 



(10) 
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Table 1. Symmetry classes of Liouville functions. 



Ry 


Rx 


Functional form 


even 


even 


P = F(x 2 ,y 2 ,z 2 ,xyz) 


even 


odd 


p = y F(x 2 ,y 2 ,z 2 ,xyz) 


odd 


xFi(x 2 , y 2 , z 2 , xyz) + zF 2 (x 2 , y 2 , z 2 , xyz) 



A natural way of expanding Liouville functions on the unit sphere is the use of spherical 
harmonics, 



p(9, <f>) = £ C lm Y{ 



'ir. 



hn 



where the angles 9 and are related to the cartesian coordinates in Eq. (^) in the usual 
way: x = sm9 cos0, etc. 

We shall use the common (but not universal) sign conventions [11, 12] 
"2Z + 1 (I -my.] 1 / 2 . . . , , ns 

(12) 



4tt (l + m)\ 



Y, 



i ^> 1 



-i) m Yr(e,0), 



where m > and P i m (cos6 l ) are the associated Legendre polynomials. The Yf- 
othonormal, so that 



a 



i in 



Yr(e,<j>) P (e,<f>)dn, 



arc 



(13) 



where dfl = sin 9 d6 dep. In particular, the total mass, namely / pdVt = yAnCoo, is constant 
for our area-preserving map. We shall henceforth ignore the trivial Coo component and 
consider only the entropy of the other ones. 

The nontrivial C\ m transform as follows: during a twist, 9 is constant, and 



4> — > 4>' = 4> + a cos9. 
We thus have p'(9, <fi + a cos 9) = p(9, (ft), or 
p'(9, (p) = p(9, — a cos9), 



(14) 



(15) 



whence 



C' 



Yj n *(6,(f>)p(9,<f)-a cos 9) dQ 



(16) 



because dfl' = dVL. Substitution of flTj] ) into (16) then gives 



(17) 



where the are components of a unitary matrix: 



u: 



T {m) 



I 



—ima cos 







dn. 



(18) 



An accurate method for computing for large j and I is described in Appendix B. 

We thus see that a twist leaves m invariant but it introduces all the I > \m\ (with 
exponentially small coefficients for large I), while a rotation mixes different m but leaves / 
invariant. This emergence of higher and higher I occurs only in the classical twist and turn 
map, which thus behaves in a quite different way from the corresponding quantum map 
[6, 7]. This is because in quantum theory / has the meaning of total angular momentum, 
and the latter is a constant of the motion under a twist (which is generated by J z 2 ), while 
the classical I has no such dynamical meaning and therefore need not be conserved. 

Note that a pure twist is a regular motion: all the orbits are closed circles around 
the z-axis. We examined the growth of S, as a function of a continuous twist angle a, 
for two maps starting with p proportional to x (odd symmetry class) or to y (even-odd 
symmetry). These are represented by initial states with 



respectively, and all other C\ m = 0. Since different values of m do not mix during a twist, 
the D-entropy is the same in both cases. Figure 3 shows the result: the growth of e s is 
roughly linear, as expected. 




(19) 



S 



We now turn our attention to the rotations of the unit sphere. With spherical harmonics 
used as a basis, the matrix representation of a 90° rotation is well known [13, 14]: the 
index I is not affected, and for each I the {21 + 1) components indexed by m undergo a 
unitary transformation. Appendix B explains how to construct accurately these unitary 
matrices (we proceeded up to / = 500). 

We checked the accuracy of our numerical calculations by verifying that unitarity held 
at each step, with an error less than 10~ 8 . To achieve this result, we had to use a range of 
values of I that increased by more than a factor 3 at each step. We thus had, at the fifth 
step, components C\ m with / up to 500. This implied that the rotation matrices had all 
possible odd orders up to 1001, and twist matrices had all orders up to 500. The next step 
would have exceeded the capacity of our computer (or entailed a severe loss of accuracy). 

Figure 4 shows how S grows with the number of steps. In the case a = 3, we considered 
two different initial Liouville densities, given by Eq. (19). These functions, which belong 
to different symmetry classes, turn out to have roughly linear rates of growth of their 
D-entropies, as expected, but, surprisingly, these rates are manifestly different from each 
other. We must therefore conclude that the growth of the D-entropy, contrary to the 
Kolmogorov-Sinai entropy [8, 9], is not related in a simple way to the Lyapunov exponents 
of individual orbits, because generic aperiodic orbits have no symmetry. 

We also performed similar calculations for some lower values of a. For a = 2, the twist 
and turn map is mostly regular, but there still are small chaotic regions [6]. These chaotic 
regions become almost invisible for a = 1.5 or less. We thus expected to find a growth of 
S starting about logarithmically, as for regular systems, and gradually becoming linear, 
with a small slope, due to the presence of small chaotic regions. However it turned out, 
as Fig. 4 shows, that for these low values of a the growth of S is not uniform. Rather, 
S oscillates about a slowly increasing average. It is therefore hopeless to detect a clear 
transition from the logarithmic regime to the linear one, and we did not proceed further 
in these numerical simulations. (We also found oscillations for a = 3, when we started 
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with an unsymmetrized /, for example with Cn = 1, and all other C/ m = 0.) 

IV. CONCLUDING REMARKS 

We found some results that were not unexpected, but we also had several surprises 
for which we can offer no explanation, and which may perhaps be worth further inves- 
tigation. There can be no doubt that, for a given Liouvillian, the rate of growth of the 
D-entropy depends on the symmetry class of the Liouville density. Therefore, contrary to 
the Kolmogorov-Sinai entropy [8, 9], the D-entropy is not directly related to the Lyapunov 
exponent of classical trajectories. 

Even more surprising is a quantitative comparison of Figures 3 and 4. The growth of 
S for a pure twist (which is a regular mapping) is faster than its growth for a weakly 
chaotic twist and turn map (with low a), when we compare the total twist angle in the 
first case, and the sum of discontinuous twists of the discrete map. 

The initial motivation of this work was a search for the existence of genuine quantum 
chaos (namely, phenomena that are not only a quantum simulation of classically chaotic 
properties). We have found definite clues that quantum chaos may indeed appear if, and 
only if, the quantum system is unbounded and has a continuous spectrum. 
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APPENDIX A. THE LYAPUNOV EXPONENT 



Consider infinitesimal deviations x — > x + e£, y — ► y + erj, and z — > z + e£, from a 
fiducial orbit of the map @. These deviations satisfy 

x£ + ^ + zC = 0, (20) 

so as to lie on the surface of the sphere. We can also define the tangential components of 
the vector (£, 77, C)> namely 

u = (—y£ + xr]) /yl — z 2 and v = Q/yl — z 2 , (21) 

in the "east" and "north" directions, respectively. 

Likewise, at the next step, let x' — > x' + e£', etc. Neglecting terms of order e 2 and 
higher, we obtain from Eq. (§), 

r = c, 

rj' = C, sin(a2;) + rj cos(az) — (az 1 , (22) 
(' = — £ cos(az) + 77 sin(a2;) + (ay'. 

This set of equations is the linear dynamical law for the evolution of the vector (£, 77, £)• 
We now want to know how an infinitesimal "unit" circle, namely a set of points with 
initial components u = cos a, v = sin a (where a runs from to 2tt), transforms into an 
ellipse (with the same area). The asymptotic growth of the major axis of this ellipse gives 
the Lyapunov exponent, as defined by Eq. (|7|). 

To find how this infinitesimal circle transforms, we note that, by virtue of the linearity 
of (^), if the initial components of a tangential vector are (cos a) and (sin a), these com- 
ponents become, after a number of steps, (Tn cos a + T i2 sin a) and (T 2 i cos a + T22 sin a), 
respectively, where the coefficients T pq are independent of a. It is therefore enough to 
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consider two infinitesimal tangent vectors having, initially, a = and tt/2. Their evolu- 
tion determines all the components of T pq , and the length of the semi-major axis is then 
obtained by maximizing the expression 



where a = J2 pq (T pq ) 2 . Finally, the Lyapunov exponent is given by Eq. (|7|), in the limit 
n — > oo. 



The evaluation of the matrix elements E/jj 71, in Eq. fll8|) is the procedure which consumes 
most of the time in our calcualtions. Each one of the indices, j, I, m, may run up to 500. 
There are therefore many millions of different integrals for which no analytic expression 
is known. These integrals are, for large values of their indices, rapidly varying functions 
of cos 8, and a straightforward numerical integration, with the level of accuracy that we 
wanted, would be prohibitive. 

We took advantage of the fact that a twist by a finite angle a can be generated by a 
sequence of infinitesimal twists by an angle e, for which we can replace e~ tmez by (1— imez). 
This entails no loss of precision if me < 10~ 16 , when we compute with 16 significant digits. 
The matrix elements U^ r l a \e) can be evaluated explicitly (see below), and each matrix is 
then raised to the appropriate power. For example, by taking e = 2~ k a, we merely have 
to compute (. . . {{U 2 ) 2 ) 2 . . .) 2 , k times. 

We only explicitly need 



r 2 (a) = (Tn cos a + T 12 sin a) 2 + (T 2i cos a + T 22 sin a) 2 , 



(23) 



as a function of a. The result after n steps is 



[ (T+ ( (7 2_ l) l/ 2 ]l/2 



(24) 



APPENDIX B. THE TWIST MATRIX 




(25) 
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This expression is readily evaluated by using the identity 

z Pp(z) = [(J - m + 1) Pp +1 (z) + (j + m) Pp_ 1 }/(2j + 1), (26) 
and the orthogonality relations for associated Legendre polynomials [15]. 

APPENDIX C. THE ROTATION MATRIX 

There are explicit formulas for the unitary rotation matrices R$ n [13, 14] However, 
when j is large, these formulas are not convenient for numerical calculations, because 
each matrix element is given as the tiny algebraic sum of a large number of huge terms 
with opposite signs. Only a few elements with m = ±j, or m close to ±j, can easily be 
obtained from the general formula. 

A much more efficient way of obtaining the R matrix for a 90° turn around the y-axis 
is to use its very definition, namely R^J X R = J z , or 

J X R = RJ Z . (27) 

With the standard representation, namely (J z ) mn = m5 mn and J x real, this gives 

[(j + m ) (j - m + 1)] 1/2 RS_ 1>n + [(j + m + 1) (j - m)} 1 ' 2 R$ +1>n = 2n R$ n . (28) 

Let us therefore define, for each j and n, a "vector" V m by 



2' 



(j + m)! (j — to)! 
(j + n)\ (j - n)\ 



1/2 



V m . (29) 



From the explicit formulas mentioned above [13, 14], it can be seen that all the components 
of V m are integers, and in particular V-j = 1 and V±-j = 2n. The recursion relation (|28f) 
thus becomes 

{j to + 1) V m _i - 2n V m + (j + to + 1) V m+1 = 0. (30) 
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(This equation also appeared as the last equation of ref. [7], unfortunately marred by a 
misprint. The calculations in ref. [7] were correct.) 

Even with the simple recursion formula (30), it is not trivial to obtain all the V m for 
j > 50 (when using double precision floating point arithmetics), because the recursion 
is unstable and small numerical errors grow exponentially. We therefore performed this 
calculation exactly, using integer arithmetics; and we checked that the resulting matrices 
were indeed unitary. 
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Captions of figures 

FIG. 1. Area-preserving projection of the hemisphere y > 0. The figure shows 20 000 
points belonging to a single chaotic orbit. The empty regions are filled by regular orbits, 
not shown here. 

FIG. 2. Distribution of Lyapunov exponents for 10 5 randomly chosen orbits. Each bin 
of the histogram has width 10~ 2 . 

FIG. 3. Growth of the .D-entropy for a continuous twist, as a function of the twist angle. 

FIG. 4. Growth of the D-entropy for the twist and turn map, for various values of the 
twist parameter a. 
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